Computation of the Drift Velocity of Spiral Waves using Response Functions 
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Rotating spiral waves are a form of self-organization observed in spatially extended systems of physical, 
chemical, and biological nature. In the presence of a small perturbation, the spiral wave's centre of rotation and 
fiducial phase may change over time, i.e. the spiral wave drifts. In linear approximation, the velocity of the 
drift is proportional to the convolution of the perturbation with the spiral's Response Functions (RFs), which 
are the eigenfunctions of the adjoint linearized operator corresponding to the critical eigenvalues A = 0, iioj. 
Here we demonstrate that the response functions give quantitatively accurate prediction of the drift velocities 
due to a variety of perturbations: a time dependent, periodic perturbation (inducing resonant drift); a rotational 
symmetry breaking perturbation (inducing electrophoretic drift); and a translational symmetry breaking pertur- 
bation {inhomogeneity induced drift) including drift due to a gradient, step-wise and localised inhomogeneity. 
We predict the drift velocities using the response functions in FitzHugh-Nagumo (FHN) and Barkley models, 
and compare them with the velocities obtained in direct numerical simulations. In all cases good quantitative 
agreement is demonstrated. 

PACS numbers: 02.70.-c, 05.10.-a, 82.40.Bj,82.40.Ck, 87.10.-e 



I. INTRODUCTION 

Spiral waves are types of self-organization observed in 
physical [1-3], chemical [4, 5], and biological [6-1 1] systems, 
where wave propagation is supported by a source of energy 
stored in the medium. The interest in the dynamics of spiral 
waves has significantly broadened in the last decade as the de- 
velopment of experimental techniques has permitted them to 
be observed and studied in an ever increasing number of di- 
verse systems such as magnetic films [12], liquid crystals [13], 
nonlinear optics [14, 15], novel chemical systems [16], and in 
subcellular [17], tissue [18] and population biology [19]. 

In the ideal unperturbed medium, the core of a spiral wave 
may be anywhere, depending on initial conditions. However, 
real systems are always subject to a perturbation. A typical re- 
sult of a symmetry-breaking perturbation is drift of the spiral 
waves, which has two components, temporal drift, which is 
shift of spiral wave rotation frequency, and spatial drift, that is 
slow movement of the spiral's rotation centre. Drift of spiral 
waves, particularly the spatial drift, is of great practical inter- 
est to applications. In cardiac tissue, drift of re-entry circuits 
may be caused by internal tissue inhomogeneities, or by ex- 
ternal perturbations, such as electrical stimulation. The possi- 
bility of control of arrhythmias by weak electrical stimulation 
has been a subject of intensive research for decades. 

Understandably, the drift of spiral waves was mostly stud- 
ied in the BZ reaction, which is the easiest excitable sys- 
tem for experimental study, and in the heart tissues and tis- 
sue cultures, which represents the most important application 
area. Examples of drift observed in experiments and numer- 



ical simualtions include "resonant" drift caused by (approx- 
imately) periodic modulation of medium properties through 
external forcing [20], constant uniform electric field that 
causes electrophoresis of charged ions taking part in the chem- 
ical reactions [21], a spatial gradient of medium properties 
[22-25] and pinning (anchoring, trapping) to a localized inho- 
mogeneity [26-28]. Interaction with a localized inhomogene- 
ity can be considered to be a particular case of the general 
phenomenon of vortex pinning to material defects, ranging 
from convective microvortex filaments in nanosecond laser- 
matter interaction to magnetic flux strings in the Sun's penum- 
bra [29]. 

A most intriguing property of spiral waves is that despite 
being propagating waves affecting all accessible space, they, 
or rather their cores, behave Uke point-Uke objects. 

Correspondingly, three-dimensional extensions of spiral 
waves, known as scroll waves, act as string-like objects. 
There have been several ad hoc theories of drift of spiral and 
scroll waves exploiting incidental features in selected mod- 
els, e.g. [30-33]. Our present study is based on an asymptotic 
theory applicable to any reaction-diffusion system of equa- 
tions in which a rigidly rotating spiral wave solutions ex- 
ist. The theory was first proposed for autonomous dynam- 
ics of scroll waves for the case of small curvatures and small 
twists [34, 35] and then extended to the drift of spiral waves 
in response to small perturbations [36]. In this theory, the 
particle-like behaviour of spirals and string-like behaviour of 
scrolls corresponds to an effective localization of so called re- 
sponse functions (RFs, see exact defininition later in this pa- 
per). The localization of RFs is the crucial assumption, which 
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underpins the entire analysis. Originally [37] this property 
was only a conjecture based on the phenomenology of spi- 
ral waves in experiments and numerical simulations [31, 38- 
40]. The analytical calculation of the response functions ap- 
pears to be infeasible. Numerical calculations in the Barkley 
model [41] and the Complex Ginzburg-Landau equation [42] 
have confirmed that indeed they are essentially localized in 
the vicinity of the core of the spiral. The asymptotic theory 
based on the response functions has been successfully used 
to quantitatively predict drift of spirals, for resonant drift and 
drift due to parametric inhomogeneity in the CGLE [43^5] 
and for drift in response to a uniform electric field in Barkley 
model [46]. Despite this success, so far the asymptotic the- 
ory has not become a generally used tool for the prediction of 
spiral wave drift. This is partly due to difficulties in the numer- 
ical calculation of the response functions. In our recent pub- 
lication [47] we have presented an efficient numerical method 
of calculating response functions in an arbitrary model with 
differentiable right-hand sides. The complexity of calculat- 
ing response functions with this method is similar to the com- 
plexity of calculating spiral wave solutions themselves. In the 
present paper, we describe the application of the asymptotic 
theory using the response functions for the prediction of sev- 
eral types of drift and show how it works for two of the most 
popular generic excitable models, the FitzHugh-Nagumo sys- 
tem [48-50], and the Barkley system [51]. We demonstrate 
that predictions of the asymptotic theory are in good quantita- 
tive agreement with direct numerical simulations. In addition, 
we demonstrate that the response functions are capable of pre- 
dicting nontrivial qualitative phenomena, such as attachment 
of spiral waves to stepwise inhomogeneity and orbital move- 
ment around a localized inhomogeneity. 

The structure of the paper is as follows. In Section II, we 
briefly recapitulate the asymptotic theory of the drift of spiral 
waves in response to small perturbation and present explicit 
expressions for drift parameters in terms of the spiral wave's 
response functions for several sorts of drift. In Section III, 
we describe the numerical methods used for calculating the 
response functions, for direct numerical simulations, and for 
processing of the results. The results are described in Sec- 
tion IV. We conclude the paper by discussion of the results 
and their implications in Section V. 



II. THEORY 
A. General 

We consider reaction-diffusion partial differential equa- 
tions. 



A rigidly clockwise rotating spiral wave solution to the sys- 
tem ( 1 ) has the form 



U = U{p{r - R),'d{r- R) + cut - 



(2) 



where R = {X, F)^ is the center of rotation, $ is the initial 
rotation phase, and p{f — R), ■d{f — R) are polar coordinates 
centered at R. For a steady, rigidly rotating spiral, R and $ 
are constants. The system of reference co-rotating with the 
spiral's initial phase and angular velocity uj around the spi- 
ral's center of rotation is called the system of reference of the 
spiral. In this system of reference, the polar angle is given by 
6* = 1? + w< - (f>, with ^ = and <3> = 0. In this frame, 
the spiral wave solution U(p, 6) does not depend on time and 
satisfies the equation 

f (U) - cjUe + DV^U = 0, (3) 

where the unknowns are the field U(p, 9) and the scalar lu. 
In a slightly perturbed problem 



dtu ^ f (u) + DV^u + eh, he 



kl « 1, 



(4) 



where eh(u, r, t) is some small perturbation, spiral waves 
may drift, i.e. change rotational phase and/or center loca- 
tion. Then, the center of rotation and the initial phase are 
no longer constants but become functions of time, R = R{t) 
and $ = In the co-rotating system of reference, time 

dependence will take form of a phase depending on time 

(j){t) ^LOt- 

Thus, we consider three systems of reference: 

1. laboratory, (f, t) ; 

2. co-moving, (p, i?, t), where (p, -d) = 

(^p{r — R), d{r — R)^ is the polar coordinate 
system centered at R\ 

3. co-rotating, (p, 6, cj)), where 9 ~ ■d{f— R) + (f){t) is the 

polar angle, and </> = cji — <I>(t) is the rotational phase, 
replacing time. 

We shall look for a solution to (4) in the form of a slightly 
perturbed steady spiral wave solution 

U(p,0» = U(p,0) + eg(p,0,0), 

where g G M*^, < e < 1. 
Then, assuming that 

i?, $ = 0(e), 

at leading order in e, the solution perturbation g will satisfy 
the linearized system 



atu = f(u) + DV2u, u.feM'', De 



>2, (1) 



where \i{r,t) = [ui^ . . .ui)^ is a column-vector of the 
reagent concentrations, f(u) ~ {fi, . . . fi)"^ is a column- 
vector of the reaction rates, D is the matrix of diffusion coef- 
ficients, and r G is the vector of coordinates on the plane. 



M^-/:)g = H(u,p,0,, 



(5) 



where 



H(U,p,0,0)=h(U,p,0,0)-- 



^R - dglJ^ 
dR 
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where h(U, p, 6, </>) is the perturbation h(u, r, t), considered 
in the co-rotating frame of reference. 
The linearized operator 

C = -DV^ -udg+d^fiV), (6) 

has critical (Re (A) = 0) eigenvalues 

£V(") = A„V("\ A„ = inw, n = 0, ±1, (7) 

which correspond to eigenfunctions related to equivariance of 
(1) with respect to translations and rotations, i.e. "Goldstone 
modes" (GMs) 

v(o) = -dgV{p,e), 

V(±^) = -ieT"'(9pTip-'9e)U(p,0). (8) 

In this paper we do not consider perturbations h(u, r, t) that 
depend on t other than 27r/aj-periodically (for a more general 
version of the theory free from this assumption see [36, 45]). 
Then h(U, p, 9, (p) is a 27r-periodic function in (j>, and we look 
for solutions g{p, 6, (p) to equation (5) with the same period- 
icity. A solvability condition leads to the following system of 
equations for the drift velocities, 

$ = .£^(w("),h(U,p,0,0))g + O(.2), 



and 



where R = X + iY is the complex coordinate of the in- 
stant spiral centre, the inner product (• , •) stands for the scalar 
product in functional space 



(w , v) = J w (r) v(r) d r, 

and the kernels W'") {p, 9), n = 0, ±1, are the response func- 
tions, that is the critical eigenfunctions 

= //„W("\ pn = -iriw, n = 0, ±1, (9) 

of the adjoint operator £+, 

C+ = DV^ + ^dg + (9uf (U))^ , (10) 

chosen to be biorthogonal 

yw, vw\ = <5,-fc, (11) 



to the Goldstone modes (8). 

The drift velocities can be written as (henceforth we shall 
drop the O(e^) terms) 

$ = eFo(i?,$), i? = eFi(i?,$), (12) 

where the "forces" i^o and Fi = (Re {Fi) , Im {Fi))^ are de- 
fined by 

F4i?, $) = (W(") (p, 9) , a„(p, 9; R, $)) , 

n = 0,l, (13) 



a„(p,6l;i?,$) 



2ir 



h(U,p,^,0)-^. (14) 

ZTT 



In the above formulae, the dependence on {R, $) is explicitly 
included to emphasize that the response functions depend on 
coordinates (p, 6) in the corotating frame of reference whereas 
the perturbations are typically defined in the laboratory frame 
of reference, and the two systems of references are related via 
R and $. 

Below we show how the forces (13), determining the ve- 
locity of the drifting spiral wave subject to a variety of pertur- 
bations, can be calculated using the computed response func- 
tions W("). We also compare the quantitative analytical pre- 
diction of drift velocities with the results of direct simulations. 



B. Resonant Drift 

Let us consider a spiral wave drifting due to the perturbation 

h(u, f,t) = A cos{ujt), (15) 

where A G is a constant vector In the co-rotating frame 
the perturbation (15) will be 



h = Acqs ((/) + <&) 

Substitution of (16) into (14) gives 

A 

OiQ = 0, ai = —e 

and, by (13), 

Fo = 0, Fi 
Hence the speed of the resonant drift of the spiral is 

e/w(i), A^ 



(16) 



1 



e**(W(i)(p,6'), A 





1 


R 






~ 2 



(17) 



(18) 



whereas its direction is constant and arbitrary. 



arg 



are- ( (W^^^ , A 



$:=0, (19) 



as it is determined by the inial phase of the spiral <I>, or, rather, 
by the phase difference between the spiral and the perturba- 
tion, (19) is only valid in the asymptotic sense, and a more 
accurate formulation is 



(20) 



Hence, at finite e the resonance is expected to be imprecise, 
and a typical trajectory of a resonantly drifting spiral is a circle 
of radius = \R\/\M = 0{e-^). 
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C. Electrophoretic Drift 

Here we consider an anisotropic perturbation which breaks 
rotational symmetry 



h(f) = B 



9U 

dx 



(21) 



where B S M^^^ is a constant matrix. This perturbation cor- 
responds to action of an external electric field on a chemical 
reaction where some of the species are electrically charged. 
In this case matrix B is diagonal and its nonzero elements 
represent motilities of the ions of the reaction species. The 
same sort of perturbation appears in the asymptotic dynamics 
of scroll waves [34, 35], where B = D. 

In the co-rotating system of reference, the perturbation (21) 
can be written using the Goldstone modes (8), as 



h(U,p,6i,( 



(22) 



which, by substituting into (14), gives 

Jo 27r 

(23) 

Thus, ao = 0, OLi{p,e) = -BV^^), which following (12) 
and (13) gives the velocity of the electrophoretic drift 



^= -e(w(i)(p,6'), BV(i)(p,( 



(24) 



which remains constant in time. 



D. Inhomogeneity induced Drift 

1. General 

We now consider the case when the reaction kinetics f in 
(1) depend on a parameter p, and the value of this parameter 
varies slightly in space, 

f = f (u, p) , p{r) = pfi + epi{r). (25) 

Substitution of (25) into (1) gives, to the first order in e, 

dtVi = DV2u + f(u,po) + epi(r)apf(u,po), 
with the perturbation in the laboratory frame of reference 

h(u,f,t) = dpi{u,pQ)pi{f). (26) 
Substitution of (26) into (14) gives 

a„(p,^?) =apf(U(p,0),po)e-"''A'„(p), (27) 

where 



(28) 



and pi {p, is the parameter perturbation considered in the 
co-moving frame of reference. The final equations for the drift 
velocities can then be written in the form 



27r OO 



(29) 





27T OO 

R = f J J w^^\p,e)c~'^Ki{p) pdpd9, (30) 



where for brevity we introduce 

«;(")(p,^)= \^^'''\p,e)Y dpi{p,e-po). (31) 



2. Linear Gradient 

Let pi vary linearly in a sufficiently large region containing 
the spiral tip and its subsequent drift trajectory. Specifically 
we consider pi ^ x — xq, where the x-coordinate of the tra- 
jectory remains near .tq. In the co-moving reference frame, 
the linear gradient perturbation will be 



Pi ^ X — xa + pcos(i?). 



(32) 



Substituting (32) into (28) gives 

Kn{p) = {X - XQ)5nfi + i/0((5„,l + . 

Then, by (27), (12) and (13), the velocity of the drift due to 
gradient of a model parameter will be 



$ = e(X-a-o) j J w^°^p,e) pdpde, 





27r OO 

R = \j j w'^^\p,e)c-'%^dpde. 





(33) 



An important feature of equations (33) is that the first of 
them depends on X while the second does not. The depen- 
dence on X means that the drift velocity changes during the 
drift, unless the drift proceeds precisely along the y-axis. As 
it happens, at first order in e, only the temporal drift, that is the 
correction to the frequency, shows this dependence. Namely, 
the first of equations (33) shows that the instant rotation fre- 
quency corresponds to the parameter value at the current cen- 
tre of rotation, p ~ pa + epi = po + ^ xq)- The spatial 
drift, described by the second of equations (33), does not de- 
pend on X. That means that while the drift proceeds, its speed 
and direction remain the same, at least at the asymptotic order 
considered. This is an important observation, firstly, because 
it allows us to treat linear gradient induced drift in the same 
way as the electrophoretic drift, i.e. expecting drift along a 
straight line, and secondly, that unlike electrophoretic drift, 
the assumption is inherently limited to such X that e{X — xq) 
remains sufficiently small. 
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3. Step inhomogeneity 

Here we consider a step perturbation located a.tx = Xs, 

Pi{x) = H(x- - Xs), 

where H() denotes the Heaviside unit step function. In the 
co-moving frame of reference we have 

=H(X + pcos(z?) -x,). (34) 

Substitution of (34) into (28) gives 

x,-X\ di? 



Kn^ I cos(nd) H 
Jo 



cos(i?) 



P / 27r 



We consider three intervals for 



-X 



(1) p <\xs~ X\,Xs> X. Then H {cos{d) - ) = 0, 

therefore iiTo = Ki = 0. 



(2) p <\x,- X\,Xs < X. Then H (^cos(i?) - ^^^^) = 1, 
therefore = 1, iCi = 0. 

(3) /9 > - X\. Then, for i?o = arccos (^^^^), 



H cos(i9) 



, otherwise. 



Thus, 



1 (xs-X 

Kq = —arccos 

TT \ p 



(35) 
(36) 



Substituting the above if„ for the three intervals into (13) and 
(12), we get the velocities of the drift due to a step- wise inho- 
mogeneity of a model parameter in the form 



2-K oo 



Xs- X^ 



R^-j J w^'\p,e)e-'" dl - i^-^^j pdpdO, (37) 

\x,-X\ 

2-K oo 2ir \xs-X\ 

^=^J J w^°\p, d) aiccos [^^— ^ p dp dd + ell{X -xs) j j w^°\p,0)pdpde. (38) 

\x^-X\ 



Note that both R and $ are functions of the current x- 
coordinate of the spiral with respect to the step, d = X — Xg, 
and R is an even function of this coordinate. 



4. Disk-shaped inhomogeneity 

We now consider an inhomogeneity which is unity within 
a disc of radius Ri„ centered at {xd,yd), and which is zero 
outside the disc. Thus we have 

= H {Rl -{x- xdf ~{y- Vdf) . 

Then calculations, similar to those for a stepwise inhomo- 
geneity, lead to 

K, = ^ arccos (^^!±|— ii^ , (39) 



where / and -do designate the distance and the direction from 
the current centre of the spiral to the centre of the inhomo- 
geneity, i.e. Xd ^ X + I cos "do, yd = y + I sin ?9o- 

This leads to the equations for the drift velocities in the 
form 
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27r l+B.i„ I T 

/ / 2 I 72 _ p2 \ ^ 



R = Ic'"" J J (^^!±^--^j pdpdf?, (41) 

\l-R,n[ 

In 2 2 2 

$ = y" iL,(")(p,6l)arccos (^^-ti^— pdpd6', (42) 

|i-i?,i„ 



It is straightforward to verify that if '&q ^ 0, Xd ~ Xs + Rm 
and Rin oo, that is when the disk is so large it turns into a 
half-plane at x > Xg, then expressions (39) and (40) tend to 
expressions (35) and (36) respectively, as should be expected. 
Another interesting limit is R„ — > 0, in which we get 



2ir 



de 



i?ir 



2ir 



in accordance with the case of a pointwise, (5-function inho- 
mogeneity considered in [52]. 



ni. METHODS 
A. Models 

We have considered two different kinetic models, both two- 
component, ^ = 2, with one nonzero diffusion coefficient. 



D 



1 




We designate u ~ f = (/,.?) 



for convenience. The kinetics FitzHugh-Nagumo system was 
chosen in Winfree [50] notation, 

g{u,v) = a{u + jv), 

with parameter values a — 0.3, /3 = 0.68, 7 — 0.5 as in [47]. 
The Barkley [51] kinetics is given by 

f{u,v) = c'^u(l - {v + b)/a), 

g{u,v) = u-v, 

with parameter values a = 0.7, b = 0.01 and c ~ 0.025, as 
in [46]. Note that both a and c are called e in [50] and [51] 
respectively; however we use e for the small parameter in the 
perturbation theory. 

B. Response functions computations 

For both the FitzHugh-Nagumo and the Barkley models, 
the response functions and the Goldstone modes have been 



computed using the methods described in [47]. The discretiza- 
tion is on disks of radii from pmax = 12 up to pmax = 50, 
using Ng = 64 of discretization intervals in the angular direc- 
tion and a varying number Np of discretization intervals in the 
radial direction, up to Np ~ 1280. The components of the spi- 
ral wave solution and its response functions for Barkley model 
are shown in fig. 1(a). Similar pictures for the FitzHugh- 
Nagumo model can be found in [47]. 



C. Perturbations 

We considered similar types of perturbations eh(u, r, t) in 
both FitzHugh-Nagumo and Barkley models, both for theoret- 
ical predictions based on response functions and in numerical 
simulations. Specifically, the perturbations were taken to have 
the following forms. 

a. Resonant drift 



h(u, f, t) = 



cos{ujt), 



(43) 



where ui is the angular velocity of the unperturbed spiral ob- 
tained as part of the spiral wave solution for the equation (3). 
b. Electrophoretic drift 



h(u, f, t) 



1 




du 

dx ' 



(44) 



c. Spatial parametric inhomogeneities. As set out in 
Section 11, a spatial dependence of a parameter p of the kinetic 
terms in the formp(r) = po +pi{r), \pi \ <C \po \ corresponds 
to the perturbation 



h(u,f,f) = dpf{u,pa)piir). 



(45) 



For each of the two models, we consider inhomogeneities in 
all three parameters, namely p G {a, /3, 7} for the FitzHugh- 
Nagumo model, and p G {a, b, c} for the Barkley model. The 
"linear gradient" inhomogeneity is of the form 



Pl = X - Xo, 



(46) 



where a^o is chosen to be in the middle of the computation box 
and close to the initial centre rotation of the spiral wave. 
The "stepwise" inhomogeneity is of the form 



pi = H(a; - Xs) 



1 



(47) 




0.624892 



u 



1.77696 



Re fW(i)l 



m 



Im (W(i)) 



(b) 



FIG. 1: (color online) (a) Solutions of the nonlinear problem (3) and the adjoint linearized problem (9,10), i.e. the response functions, as density 
plots. Barkley model, pmax ~ 12.8, Np = 1280. Numbers under the density plots are their amplitudes A: white of the plot corresponds to 
the value A and black corresponds to the value ~A of the designated field. Upper row: 1st components (u), lower row: 2nd components (v). 
The central areas ofW("',n = 0,1, are also shown magnified in the small corner panels, (b) Snapshot of spiral wave in the Barkley model 
(u: red colour component, v: blue colour component), drifting in a stepwise inhomogeneity of paramer c (green colour component). The thin 
white line is the trace of the tip of the spiral in the course of a few preceding rotations. Yellow circles are positions of the centres calculated as 
period-averaged positions of the tip. 



where Xg is varied and chosen with respect to the initial centre 
of rotation of the spiral wave. The — ^ term is added to make 
the perturbation symmetric (odd) about x = Xs,lo minimize 
the inhomogeneity impact on the spiral properties while near 
the step. As it can be easily seen, within the asymptotic theory, 
this term only affects the frequency of the spiral but not its 
spatial drift. 

The "disk-shape" inhomogeneity is of the form 



Pi = H(i?i, 



(48) 



where the position of the centre of the disk rd = {xd, yd)^ is 
varied and chosen with respect to the initial centre of rotation 
of the spiral wave. 



D. Drift simulations 

Simulations have been performed using forward Euler 
timestepping on uniform Cartesian grids on square domains 
with non-flux boundary conditions and five-point approxima- 
tion of the Lapacian. The space discretization step Ax has 
been varied between Ax = 0.03 and Ax = 0.1, and time dis- 
cretiation step At maintained as At = |Ax^. The tip of the 
spiral is defined as the intersections of isolines u{x, y) = 
andt;(x,y) ~ w*, and the angle of Vu at the tip with respect to 
X axis is taken as its orientation. We use {u^,,v^) = (0,0) for 
the FitzHugh-Nagumo model and (u*, w*) = (1/2, a/2 — b) 
for the Barkley model. 



E. Processing the results 

For coarse comparison, we use the trajectories of the instan- 
tanous rotation centre of the spiral wave. They are directly 
predicted by the theory. In simulations, they are calculated by 
averaging the position of the tip during full rotation periods, 
defined as the intervals when the orientation makes the full 
circle (— tt, tt], see fig. 1(b). 

For finer comparison, we fit the raw tip trajectories, i.e. we 
use theoretical predictions including the rotation of the spi- 
ral. That is, if the theory predicts a trajectory of the centre 
as i? = R{t; A,B,...) G C (a circle for resonant drift and 
a straight line for electrophoretic or linear gradient inhomo- 
geneity drifts) depending on parameters A,B,... to be iden- 
tified, then the trajectory of the tip is assumed in the form 
i?tip(t) = R{t; A,B,...)+ i?,o,ee*("*+®") where R,„,, e M 
is the tip rotation radius, w S M is the spiral rotation frequency 
and Oo G K is the initial phase. The parameters i?core, w and 
Qq are added to the list A,B,... of the fitting parameters. 



IV. RESULTS 

A. Simple drifts 

Fig. 2 shows a comparison between the theoretical predic- 
tions for the simple drifts and the results of direct numerical 
simulations of various perturbation amplitudes e. The simple 
drifts include the resonant drift, the electrophoretic drift and 
the drift in the linear parametric gradient with respect to one 
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FIG. 2: (color online) Drift speeds as functions of corresponding perturbation amplitudes. Top row: FitzHugh-Nagumo model. Bottom row: 
Barkley model. First column: resonant drift. Second column: electrophoretic drift. Third column: drift in linear gradient inhomogeneity, 
namely (c) with respect to parameter /3 in the FHN model, and (f) with respect to parameter a in the Barkley model. In the second and the 
third columns, the symbols represent simulations and the lines represent theoretical predictions. Numerical parameters: (a) Aa:: — Ap — 0.1, 
Pmax = 50, (b) Ax = Ap ^ 0.1, pmax = 25, (c) Aa; = 0.08, Ap = 0.02, pn^ax = 25, (d) Ax ^ Ap = 0.05, p,nax = 25, (e) 
Ax = Ap = 0.02, p„,ax = 12.5, (f) Aa: = Ap = 0.06, p^ax = 24. 



arbitrarily selected parameter. 

For the resonant drift, the motion equations given by (18), 
(19) and (20) can be summarized, in terms of complex coor- 
dinate i? = X + iF, as 

di? i* d<i> 

-^=e p, - = q, (49) 

where P — ^ |e(W'^-' , A)| is predicted by the theory at 
leading order, and q = 0{e^) is not, and we only know its 
expected asymptotic order. The theoretical trajectory is a cir- 
cle of radius p/q, and the spiral drifts along it with the speed 
p. In the simulations, we determined both the radius and the 
speed by fitting. The speed is used for comparison and the 
radius is ignored. 

For the other two types, electrophoretic drift and linear gra- 
dient inhomogeneity drift, the theory predicts drift at a straight 
line, according to (24) and (33) respectively. In these cases, 
we measure and compare the x and y components of the drift 
velocities separately. 

For numerical comparison in the case of Unear gradient in- 
homogeneity, we chose a pieces of trajectories not too far from 
X — xq, selected empirically to achieve a satisfactory quality 
of fitting. 

A common feature of all graphs is that at small enough e, 
there is a good agreement between theory and simulations. As 
expected, differences appears for larger e with the disagree- 
ment occurring sooner (for smaller values of drift speed) for 



the linear gradient inhomogeneity drift. This is related to an 
extra factor specific to the inhomogeneity-induced drift: the 
properties of the medium where they matter, i.e. around the 
core of the spiral, changes as the spiral drifts. Since we re- 
quire a certain number of full rotations of the spiral for fit- 
ting, faster drift meant longer displacement along the x axis 
and more significant change of the spiral properties along that 
way, which in turn affects the accuracy of the fitting. 



B. Numerical convergence 

Fig. 3 illustrates numerical convergence of results with dis- 
cretization parameters. We consider the simple drift cases 
and focus on forces, defined as the drift speed/velocity per 
unit perturbation amplitude e. The discretization parameter 
that primarily dictates the accuracy of solutions is a spatial 
discretization step: Aa; in the simulations and the radial dis- 
cretization step Ap in the response functions calculations. 

In simulations, the forces are determined for values of e 
well within the linear range as determined in Fig. 2. These are 
calculated for different values of the space discretization step 
Ax, where the time discretization changed simultaneously so 
that the ratio At/{Ax)'^ remaines constant. 

In theoretical predictions, the forces are given by the val- 
ues of the corresponding integrals of response functions as 
described by Section II, and we have calculated the response 
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FIG. 3: (color online) Numerical convergence of drift forces. Top row: FHN. Bottom row: Barkley. First column: resonant drift. Second and 
third columns: x and y components of the electrophoretic drift. The forces were measured for (a) e = 10"^, (b,c) e = 10"^ (d) e = 5 ■ 10"^, 
(e,f) e = 10"^ 



functions and the corresponding integrals with various values 
of the radius discretization steps Ap. 

Our discretization in both the theoretical and stimulations 
cases is second order in Aa; and in Ap, so one would expect 
to see linear dependence of the drift forces on the squares of 
these discretization steps, (Ax)^ and (Ap)^, at least for the 
values of these steps small enough. This is indeed what is 
observed. 

We have gone further and extrapolated the calculated the- 
oretical and simulation values of forces to zero Ap and Ax 
respectively, based on the expected numerical convergence 
properties. Such extrapolation gives the values of the forces 
which differ from the exact value only due to other, smaller 
discretization eiTors, which are: angular discretization and re- 
striction to the finite domain in the theoretical predictions, and 
second-order corrections in e and the boundary effects in the 
simulations. Comparison of such extrapolated data shows a 
very good agreeement between theory and direct numerical 
simulations (DNS) which is illustrated in Table 1. Note that 
the values for fig. 3(e) and fig. 3(f) are also in good agreement 
with the results of [46]. 

For the extrapolation, we fitted the numerical data with the 
expected numerical convergence dependencies, which were 
different for theoretical calculations and for the simulations. 
In simulations, the central difference approximation of the 
Laplacian means that the next term after (Ax)^ is (Ax)^. 
The expected error due to time derivative discretization is a 
power series in At ~ (Ax)^, hence the next term there after 
(Ax)^ is again (Ax)^. The situation is different in the re- 
sponse functions calculations as there is no symmetry in the 



approximation of p-derivatives, therefore we expect that in the 
theoretical convergence, the next term after (Ap)^ is (Ap)^. 
We note, however, that approximation of both theoretical and 
simulation data with similar dependencies, be it with a cubic 
or a quartic third term, gave very similar results. 



C. Drift near stepwise inliomogeneity 

The theoretical predictions for stepwise inhomogeneity, 
(45,47) and disk-shaped inhomogeneity considered next, are 
more complicated than the simple forms of drift considered 
up to this point. Because now the medium is inhomogeneous 
in the presence of the perturbation, the velocity depends on 
the instant position of the spiral center and as a result the spi- 
ral trajectories can be quite complex. Qualitative comparisons 
between theory and simulations can be made in the general 
case, but for detailed quantitative comparison we focus on 
the cases where the theory predicts simple attractors, e.g. a 
straightforward drift along the step for the stepwise inhomo- 
geneity. 

Equations (37) give a system of two first-order autonomous 
differential equations for X = Re (R) and Y ~ Im (R), 

X = eF,(X-Xs), 

Y = eFyiX~Xs), (50) 
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FIG. 4: (color online) Drift in stepwise inhomogeneity (45,47). First row: theoretical predictions for the drift forces components as functions 
of the distance to the steps, d = X — Xs, in parameters (a) a, (b) P and (c) 7, in the FitzHugh-Nagumo model. Second row: same for 
Barkley model, steps in parameters (d) a, (e) b and (f) c. Third row: comparison of theoretical predictions with DNS. (g) A phase portrait 
of the drift in the FitzHugh-Nagumo model, in theory, (37), and DNS, (4,45,47), with a step inhomogeneity of parameter a (corresponds to 
panel (a)), at e = 10"^. Shown are the theoretical vector field (black arrows; the lengths are nonlinearly scaled for visualization), a selection 
of theoretical trajectories (red filled circles) and a selection of numerical trajectories (blue open circles) of the centres of the spiral waves. 
Trajectories are arbitrarily shifted in the vertical direction for visual convenience. Dashed-dotted vertical lines correspond to the root of the 
theoretical horizontal component of the speed, and the location of the step X — Xs — 0. (h) Speed of the established vertical movement along 
the stepwise inhomogeneity as in panel (g), as a function of inhomogeneity strength, (i) A phase portrait of the drift in the Barkley model with 
a step inhomogeneity of parameter c (corresponds to panel (f)), at e = 3 ■ 10"*. Notation is the same as in panel (g). 
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TABLE I: Fitting of numerical convergence of theoretical and simulation data 



where 

27r oo 
\d\ 

(51) 

27r oo 

= i J J Im(ii;(i'(p,6l)c-"') p'^ ~ d? ApdO. 

\d\ 

(52) 

The right-hand sides of system (50) depend only on X but 
not on Y, that is, the system is symmetric with respect to 
translations along the Y axis. For this reason, the roots 
: Fx{d^.) = provide invariant straight lines along the 
y-axis. An invariant line {{xg + d*, £ R} will be sta- 

ble if eF'^{d^) < and unstable if eF^{d^) > 0. Note that the 
stability of invariant lines reverses with a change of sign of e 
and also that (d) is an even function. Hence if e ^ and 
7^ then either {xg + d^,,Y} or {xg — d*, F} will be an 
attracting invariant set. 

Fig. 4(a-f) show the theoretical predictions for the drift 
forces, i.e. velocity components per unit perturbation magni- 
tude e, Fx{d) and Fy{d), on the distance d = X — Xg from 
the instant spiral centre to the step. This is done for both 
FitzHugh-Nagumo and Barkley models, for steps in each of 
the three parameters in these models. The roots of Fx{d) are 
specially indicated. One can see from the given six examples, 
that existence of roots of F.^i) is quite a typical, albeit not a 
universal, event. 

The qualitative predictions of the theory about a stable in- 
variant line are illustrated by fig. 4(g) where we present results 
of numerical integration of the ODE system (50) and the re- 
sults of direct numerical simulation of the full system. In the 
example shown, the positive root k, 2.644 of F^ is stable 
and the negative root — c?» is unstable. Hence the theoretical 
prediction for different initial conditions are: 

for X{Q) > Xg — (i» and not too big, the spiral wave will 
approach the line X = Xg+d^, and drift vertically along 
it with the speed eFy{d^) w 0.8468e; 

for X{0) < df, the spiral wave will drift to the left with ever 
decreasing speed, until its drift is no longer detectable; 

for big I X (0) I , the drift will not be detectable from the outset. 



As seen in fig. 4(g) this is indeed what is observed, both 
for the theoretical and for the DNS trajectories, and the visual 
similarity between theoretical and DNS trajectories is an il- 
lustration of the validity of the qualitative predictions of the 
theory. 

Since the generic drift is non-stationary, a quantitative com- 
parison for typical trajectories is difficult. However, the drift 
along the stable manifold X = Xg +d^ is stationary with ver- 
tical velocity given by eFy{d^) so a comparison is easily made 
using the same methods as in the case of "simple" drifts con- 
sidered in the previous subsections. The results are illustrated 
in fig. 4(h). As expected, we see good agreement between the 
theory and the DNS for small e. 

The phenomenological predictions are different for the case 
when Fx (d) has no roots, or when its roots are so large that 
\Fy{d^.) \ is so small that the drift cannot be detected in simu- 
lations. In such cases, the theoretical predictions for different 
initial conditions are: 

for |X(0)| not too large, the spiral wave will move with vary- 
ing vertical velocity component but always in the same 
horizontal direction (to the right if eFx (0) > 0), even- 
tually with ever decreasing speed, until its drift is no 
longer detectable; 

for |X(0)| too large, the drift will not be detectable from the 
outset. 

This prediction is confirmed by simulations, as illustrated in 
fig. 4(i), where we have chosen the case of inhomogeneity in 
parameter c in Barkley model, for which the smallest positive 
root is (i» « 2.867, which gives i^y((i*) « 0.1632. This value 
should be compared to F^(0) w 48.42 and ^^(0) w 119.4. 
Note also that to get the drift velocities. Fx and Fy should 
be multiplied by e which should be much smaller than cq = 
0.025. So when the spiral is further than |X — ~ 2 from 
the step, the drift is very slow and hardly noticeable, even 
though according to the theory, there should be stable vertical 
drift around X ^ Xg + d^,, which is too slow to be observed 
in normal simulations. 



D. Drift near disk-shape localized inhomogeneity 

The theoretical predictions for the disk-shaped inhomo- 
geneity (45,48), are more complicated but also more interest- 
ing. The theoretical spiral motion equation (41) has a rota- 
tional, rather than the translational symmetry of the stepwise 
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FIG. 5: (color online) Drift around disk-shape inhomogeneity (45,48) of radius Rin — 0.56. First row: theoretical predictions for the drift 
speed components as functions of the distance to the disk centre, / = ((X — Xd)^ + (Y — yd)^)^^^ , for inhomoegenity in parameters (a) a, 

(b) /3 and (c) 7, in the FitzHugh-Nagumo model. Second row: same for Barkley model, steps in parameters (d) a, (e) b and (f) c. Third row: 
comparison of theoretical predictions with DNS. (g) Angular speed of the established orbital movement around the inhomogeneity site as on 
panel (i), as a function of inhomogeneity strength, (h) A phase portrait of the drift in the Barkley model in theory, (41), and DNS, (4,45,48), 
with disk-shape inhomogeneity (green) of parameter b (corresponds to panel (e)), at e = 10^^. Shown are the theoretical vector field (black 
arrows; the lengths are nonlinearly scaled for visualization), a selection of theoretical trajectories (red filled circles) and a selection of numerical 
trajectories (blue open circles) of the centres of the spiral wave. Dash-dotted circles coreespond to the roots of the theoretical radial component 
of the drift force, (i) A phase portrait of the drift in the FitzHugh-Nagumo model with inhomogeneity of parameter 7 (corresponds to panel 

(c) ), e — 0.3. Notation is the same as in panel (h). 
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inhomogeneity. In polar coordinates {l,'&o) centered at the 
center of the inhomogeneity, so that R = Xd + iyd + /e"'", 
equation (41) can be rewritten in the form 



I = -eFril), 
Wo - eFail), 



(53) 



where and Fa are the radial and azimuthal components of 
the drift force, given by 

|i--R,i„ 



Pdpde, (54) 



27T l + R,„ 

Fa^ J J lm(w^'\p,0)c~'' 

\l-Ri„\ 



2lp 



pdpde. (55) 
(56) 



The minus sign in the first equation of (53) comes from the 
fact that in (41), the origin was placed at the instant rotation 
centre of the spiral, and the position of inhomogeneity is de- 
termined with respect to it, where as now we do the other 
way round: the origin is at the centre of inhomogeneity and 
the current position of the spiral rotation centre is determined 
with respect to it. 

In the system (53), the axial symmetry is manifested by the 
fact that the right-hand sides of (53) depend on I but not on 
•Oq, and the equation for / is a closed one. Hence roots of 
Fr{l*) = represent invariant sets, which in this case are 
circular orbits. The movement along those orbits will have a 
linear speed eFa{l*) and angular speed O = eFa{U)/U. The 
stability of these orbits is determined by the sign of eF^l^,): 
stable for positive and unstable for negative. Unlike the case 
of the stepwise inhomogeneity, now we do not have any mir- 
ror symmetries, as only positive I make sense, therefore for a 
given root U a stable circular orbit is guaranteed only for one 
sign of e but not the other. 

Fig. 4(a-f) show the theoretical predictions for the drift 
forces Fr{l) and Fa{l)- This is done for both FitzHugh- 
Nagumo and Barkley models, for inhomogeneities in each of 
the three parameters in these models. The roots of Fr{l) are 
specially indicated. One can see from the six given examples 
that existence of roots of Fr{) is rather common and often 
there is more than one root, Ij, such that = Zo < < '2 < 

The qualitative predictions of the theory are illustrated in 
fig. 5(h,i) where we present results of numerical integration 
of the ODE system (53) and the results of direct numeri- 
cal simulations of the full system. In the example shown in 
fig. 5(h), the predictions are given by fig. 5(e), which say 
that the smallest orbit has radius li w 3.724, with the or- 



bital speed Fa{li) ~ 0.003938, which is rather small com- 
pared to max(|^V(OI ~ 3.458 and max(|i^a(OI ~ 8.534 and 
hardly observable in numerical simulations. Hence in this 
case, the radial component of the drift speed Fr{l) is effec- 
tively constant-sign, and for negative e, one should observe 
repulsion of the spiral wave from the inhomogeneity until it is 
sufficiently far from it, / 3, to stop feeling it, and for pos- 
itive e, the spiral wave will be attracted towards the centre of 
inhomogeneity from any initial position I ~ 3. This is indeed 
what is observed in simulations shown in fig. 5(h) where the 
case of e > is shown, and the centre of the inhomogeneity, 
I = Iq, is attracting for the spiral wave. 

In the example shown in fig. 5(i), the inhomogeneity cen- 
tre /q = is repelling. Instead, the first orbit of radius li « 
1.7722 is attracting. The perturbation amplitude e ~ 0.3 in 
this case is quite large and comparable with the value 70 — 0.5 
of the perturbed parameter itself. We see that although the nu- 
merical correspondence between theory and DNS in this case 
is not very good (note the distances between the open circles 
and between the filled circles), the qualitative prediction of 
orbital movement remains impeccable. As expected, the nu- 
merical correspondence becomes good for smaller values of 
e, see fig. 5(g). 



V. DISCUSSION 

We have considered symmetry breaking perturbations of 
three different kinds; time-translation symmetry breaking, 
that is homogeneous in space and periodic in time ("reso- 
nant drift"); rotational symmetry breaking through differen- 
tial advective terms ("electrophoretic drift"); and spatial trans- 
lation symmetry breaking through space-dependent inhomo- 
geneities ("inhomogeneity induced drift"). The latter type in- 
cludes three sub-cases cases: a linear parametric gradient, a 
stepwise parameter between to half plane, and a parameter in- 
homogeneity localized within a disk. 

a. Quantitative: drift velocity. We have demonstrated 
that asymptotic theory gives accurate predictions for spiral 
drift: in some cases the discrepancy between the theory and 
the direct simulations was as low as 0.02%. The discrepancy 
is affected by the numerical discretization parameters, both for 
the direct simulations and for response function computations, 
and by the magnitude of the perturbation. 

b. Qualitative: attachment and orbiting. In the more 
complicated cases of spatial inhomogeneity, the response 
functions allow us to predict qualitatively different regimes of 
spiral motion, which we have been able to confirm by direct 
simulations. 

In the presence of a stepwise inhomogeneity, the centre of 
spiral wave rotation may either be attracted to one side of the 
step where it gradually "freezes", or it may get attached to the 
step and drift along it with the constant velocity. In the lat- 
ter case, the speed of the drift is proportional to the inhomo- 
geneity strength, whereas the distance at which the attachment 
happens, does not not depend on the inhomogeneity strength 
at leading order If the sign of inhomogeneity is inverted, the 
attachment occurs on the opposite side of the step and pro- 
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ceeds in the opposite direction. 

In disk-shape inhomogeneity, the situation is somewhat 
similar but more interesting. The spiral wave may be attracted 
towards the centre of the disk, or repelled from it. It may 
also be attracted to or repelled from one or more circular or- 
bits. The drift velocity along the orbits is proportional to the 
strength of the inhomogeneity, whereas the radii of orbits do 
not depend on it at leading oder. The repulsion changes to 
attraction and vice versa, with the change of the sign of the 
inhomogeneity. 

c. Prevalence of attachment and orbiting The possibili- 
ties of attachment to the step inhomogeneity and orbital move- 
ment for the disk-shape inhomogeneity are both related to the 
change of sign of the integrals of the translational response 
functions, which in turn are possible due to changes of sign 
of the components of those response functions. Not supris- 
ingly, there is a certain correlation between these phenomena. 
The graphs fig. 5(a-f) may be viewed as deformed versions 
of the corresponding graphs fig. 4(a-f). Respectively, posi- 
tive roots of Fx{d) in fig. 4(a,d,e,f) have corresponding roots 
of Fr{l) in fig. 5(a,d,e,f). However, the integrals in equations 
(51) and (54) are only similar but not identical, and the above 
correspondence between the roots is not absolute: the roots 
of Fr{l) in fig. 5(b,c) and the smaller roots of this function in 
fig. 5(a,f) have no correspondences in fig. 4. Overall, based on 
results considered, orbital motion around a localized inhomo- 
geneity seems to be more prevalent than attachment to a step- 
wise inhomogeneity. Moreover, the typical situation seems to 
be that there are multiple stationary orbits around a disk in- 
homogeneity. We have already discussed this situation in our 
recent preliminary short communication [52] where we have 
also illustrated how for the initial conditions between two sta- 
tionary orbits, the spiral wave launched into one orbit or the 
other depending on the sign of the inhomogeneity. 

The possibility of orbital drift, related to a change of sign 
of an equivalent to the function Fr{l), has been discussed at a 
speculative level in [53]. The sign change of translational re- 
sponse functions was observed in oscillatory media described 
by CGLE [54, 55]. The examples we consider here suggest 
that this theoretical possibility is in fact quite often realized 
in excitable media, and even multiple orbits are quite typi- 
cal. Theoretical reasons for this prevalence are not clear at 
present. As stated in [52], the prevalence of multiple orbits 
may be understood in terms of asymptotic theories involving 
further small parameters. So, the version of kinematic theory 
of spiral waves suggested in [56] produces an equivalent of 
response functions, which is not only quickly decaying, but 
also periodically changing sign at large radii, with an asymp- 
totic period equal to the quarter of the asymptotic wavelength 
of the spiral wave. Other variants of the kinematic theory, 
e.g. [57, 58] did not, to our knowledge, reveal any such fea- 
tures on a theoretical level. However, numerical simulations 
of kinematic equations presented in [57] showed attachment 
of spiral waves to non-flux boundaries, which in a sense is 
similar to attachment to stepwise inhomogeneity. On a phe- 
nomenogical level, such attachment is, of course well known 
since the earliest simulations of excitable media, e.g. [38]. 



d. Orbiting drift vs other spiral wave dynamics. Prop- 
erties of the orbital drift resemble properties of resonant drift 
when the stimulation frequency is not fixed as in the exam- 
ples above, but is controlled by feed-back [59]. In that case, 
the dynamics of the spiral wave is controlled by a closed au- 
tonomous system of two differential equations for the instant 
centre of rotation of the spiral, like (50) or (53). In particular, 
depending on the detail of the feedback, this planar system 
may have limit cycle attractors, dubbed "resonant attractors" 
in [60], which may have circular shape if the system with the 
feedback has an axial symmetry. Apart from this being a com- 
pletely different type of drift, we also comment that the second 
order ODE system is an approximation subject to the assump- 
tion that the feedback is instant, and in the situations when the 
delay in the feedback is significant due to the system size and 
large distance between spiral core and feedback electrode, the 
behaviour becomes more comphcated. 

For some combination of parameters, the trajectory of an 
orbiting spiral may also resemble meandering and may be 
taken for this in simulations or experiments. So, it is possi- 
ble that orbital movement was actually observed by Zou et al. 
[61, p. 802] where they reported spiral "meandering" around a 
"partially excitable defect"; although it is difficult to be cer- 
tain as no details are given. The difference is that spiral mean- 
dering, in the proper sense, is due to internal instabilities of a 
spiral wave, whereas orbital motion is due to inhomogeneity. 
E.g. in orbiting, the "meandering pattern" determined by H/lu 
will change depending on the inhomogeneity strength. 

The phenomenon of "pinning" of spiral waves to localized 
inhomogeneities has important practical implications for the 
problem of low voltage defibrillation [62-65]. In terms of 
spiral wave dynamics this is usually understood as attraction 
of the spiral centre towards the inhomogeneity locus. Practi- 
cally interesting cases of pinning are usually associated with 
inexcitable obstacles, which are not small perturbations and 
therefore not amenable to the asymptotic theory considered in 
this study. However, the possibility of orbital motion around 
a weak inhomogeneity suggests that a similar phenomenon 
may be observed in strong inhomogeneities as well. This of- 
fers an unexpected aspect on the problem of pinning. Instead 
of a simplistic "binary" viewpoint, that a local inhomogene- 
ity can either be attractive, which is the case of pinning, or 
repelling, which is the case of unpinning, there is actually a 
third possibility, which can in fact be more prevalent than the 
first two, namely, that at some initial conditions the spiral may 
orbit around one of a number of circular orbits, regardless of 
whether or not it is attracted to the center, which can be con- 
sidered just as one of the orbits that happens to have radius 
zero. That is, there is more than one way that a spiral may be 
bound to inhomogeneities. 

e. Conclusion. We have demonstrated that the asymp- 
totic theory of spiral wave drift in response to small perturba- 
tions, presented in [36, 45], works well for excitable media, 
described by FitzHugh-Nagumo and Barkley kinetics mod- 
els and gives accurate quantitative prediction of the drift for a 
wide selection of perturbations. 

The key objects of the asymptotic theory are the response 
functions, i.e. the critical eigenfunctions of the adjoint lin- 



15 



earized operator. The RFs have been found to be localized in 
most models where they have been calculated; however there 
are counterexamples demonstrating that surprises are possi- 
ble [66]. Physical intuition tells that for the response function 
to be localized, the spiral wave should be indifferent to dis- 
tant perturbations, which will be the case if the core of the 
spiral is a "source" rather than a "sink" in the sense of the 
flow of causality, for example as defined by the group veloc- 
ity. Indeed, this localization property has been proven for one- 
dimensional analogues of spiral waves [67] and there is hope 
that this result can be extended to spiral and scroll waves. 

The effective spatial localization of the RFs on the mathe- 
matical level guarantees convergence of the integrals involved 
in asymptotic theory, and on the physical level explains why 
wave-like objects like spiral and scroll waves, while stretch- 
ing to infinity and synchronizing the whole medium, behave 



respectively as particle-like and string-like localized objects. 
This macroscopic dissipative wave-particle duality of the spi- 
ral waves has been previously demonstrated for the Complex 
Ginzburg-Landau equation [45] which is the archetypical os- 
cillatory media model. Here we confirmed it for the most pop- 
ular excitable media models important for many applications. 
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